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Abstract. We consider the problem of exact reconstruction of univariate functions witli jump discontinu- 
ities at unknown positions from their moments. These functions are assumed to satisfy an a priori unknown 
' linear homogeneous differential equation with polynomial coefficients on each continuity interval. There- 

' fore, they may be specified by a finite amount of information. This reconstruction problem has practical 

, importance in Signal Processing and other applications. 

1 It is somewhat of a "folklore" that the sequence of the moments of such "piecewise D-finite" functions 

satisfies a linear recurrence relation of bounded order and degree. We derive this recurrence relation ex- 
plicitly. It turns out that the coefficients of the differential operator which annihilates every piece of the 
function, as well as the locations of the discontinuities, appear in this recurrence in a precisely controlled 
manner. This leads to the formulation of a generic algorithm for reconstructing a piecewise D-finite function 
from its moments. We investigate the conditions for solvability of the resulting linear systems in the general 
case, as well as analyze a few particular examples. We provide results of numerical simulations for several 
, types of signals, which test the sensitivity of the proposed algorithm to noise. 
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1. Introduction 

Consider the problem of reconstructing an unknown function 5 : [a, 6] ~> M from some finite number of its 
power moments 

i (1-1) mkig)^ / x*'g{x)dx k = 0,l,...,M 

^ ■ . . . 

This formulation is a "prototype" for various problems in Signal Processing, Statistics, Computer Tomogra- 

■ phy and other areas (see [T1[T71[57] and references there). In all practical apphcations, it is assumed that g 
0^ , belongs to some a priori known class, and it can be faithfully specified by a finite number of parameters in 

■ that class. For example, smooth signals may be represented as elements of some finite-dimensional Hilbert 
space and then the reconstruction problem is analyzed in the classical framework of linear approximation - 

, J see [32I [1] for thorough expositions. 

rS 

^ . In recent years, novel algebraic techniques for moment inversion are being developed. One reason for their 

appearance is the unsatisfactory performance of the classical approximation methods when applied to irregu- 
lar data. The famous "Gibbs effect" due to a jump discontinuity probably provides the best-known example 
of such undesired behaviour. In [14 it is shown that only nonlinear methods have a chance to achieve the 
same order of approximation for such discontinuous signals as conventional (linear) methods do for smooth 
signals. The various nonlinear methods developed recently include the framework of signals with finite rate of 
innovation ([MIISIIS]), Pade-based methods ([TUlllSllS]), and other algebraic schemes ([TTllllll])- Methods 
for reconstructing planar shapes from complex moments are developed in '18J [17] . In [21] an approach based 
on Cauchy-type integrals is presented. In [27 ^ . several of the above methods are reviewed in detail, while 
emphasizing the mathematical similarity between the corresponding inversion systems. 

The present paper develops a general method for explicit inversion of the moment transform for functions 
which are piecewise solutions of linear ODEs with polynomial coefficients (precise definitions follow). In 
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particular, the method allows to locate the discontinuities of the function by purely algebraic means. The 
algebra involved is in the spirit of holonomic combinatorics ([35]). We believe that the tools developed 
in this paper may shed a new light on the similar structure of the algebraic equations which appear in 

[31 El HSl [13 117]. 



1.1. Overview of main results. Let £> denote an arbitrary linear differential operator with polynomial 
coefficients 



N 



(1.2a) T> = Y,PA^)d' 

3=0 
kj 

(1.2b) Pj{^) = ''^^0'i,j^^ O^i.j G 



i=0 

where d is the differentiation operator with respect to x: d^f =^ ^/ = f^^^K^) and 9" = I, the identity oper- 
ator. If some function g satisfies D g = 0, we say that J) annihilates g and also write g G A/j) =^ {/ : S) / = 0}. 
These functions are called "Z?-finite" ( "differentiably- finite" , [30|). 

The class A-^ of "piecewise 13- finite" functions is defined in the following way: for n 0, 1, . . . , /C, let 

An [^n, ^n+i] be a partition of [a, b] such that — oo < a = < Ci < ■ ■ • < — b < +oo (the case K- — 
corresponds to functions consisting of a single "piece"). We say that g £ A^i if there exist operators S)„ such 
that on each "continuity interval" A„, the n-th piece of g equals to, say, gn{x), such that 'Sngn = 0. In the 
most general setting, the annihilating operators D„ may be pairwise different. In this paper we explicitly 
treat functions for which Dn = 2D for all < n < K,. We write g G A'^ in this particular case. 

We shall always assume that the leading coefficient of D„ does not vanish at any point of A„. 

We are interested in solving the following 

Operator-Based Moment Reconstruction Problem. Given the sequence of the moments (|l.ip of an 
unknown function g G A'q and the constants /C, iV, {kj}, a, b: 

1) Find the operator D (in the general case, the operators S)„) - that is, determine the coefficients Uij 
(respectively aij^n) as in ()1.2I) such that S)„ 5„ = for all < n < /C. 

2) Find the "jump points" {^n}n=i- 

3) Let {ui}^^ be a basis for the linear space TV©. Find constants ai^n such that gn{x) ~ '^iLiCti,nUi{x). 
That is, determine each concrete solution of 'Dn9n = on A„. 

Example 1.1. Let g{x) = ae^^ be our unknown function on [0, 1]. The parameters of the problem are: 
/C = 0, D = 9 — /3I, [a, 6] = [0, 1] and the unknowns are: a, (3. Direct calculation of the moments yields 

J a{e^ - 1) = f3mo 
1 ae^ = /3mi + niQ 



(1.3) 



Denoting 2ii =^ A, ([O]) amounts to: 



(1.4) 



' \^A 



This system does not have an explicit analytic solution. It can be solved numerically (for example, using 
Newton's method) and a unique solution exists provided A G (0, 1) (the graph of the function y = ^-e-^^ ~ x 
is monotone in whole M and < y{x) < 1). ■ 



The above solution is unsatisfactory for several reasons. First, it is not general enough. Second, the solution 
is available only as an approximation and not in a closed form. However, there is an important positive 
feature: the minimal possible number of measurements is used. Using our method. Example 11.11 will be 
solved later in a more convenient and general way - see Examples 12.21 and 13.11 

Our method is based on the following results which we prove below (Sections 12. li I2.2 land l3.1|) . These results 
establish explicit relations between the known and the unknown parameters of the reconstruction problem. 

Theorem (|2.9l 13. ip . Let JC = and S) g = 0. Then the moment sequence of g satisfies a linear recurrence 
relation with coefficients linear in atj. Consequently, the vector a= (oij) satisfies a linear homogeneous 
system Ha. = where the entries of H are certain linear combinations of the moments of g whose coefficients 
depend only on N and the endpoints a, h. 

Theorem (|2.12p . Let /C > and let "S) annihilate every piece of g & A^. Then the operator given by 
={]Ai=i{x ~ it)^ annihilates (7 as a distribution. Consequently, conclusions of Theorems \ 2.9\ and 
\3.1\ are true with S) replaced by S) . 

Proposition (j3.3p . Let g G with operator £) annihilating every piece g„. Let {ui]f^i be a basis for the 
linear space Af's and gn{x) — '^f^i C(i,nUi{x) . Then the vector a of the coefficients ai^n satisfies a linear 
system Ca — m where the matrix C contains the moments of Ui and the vector m contains the moments of 
9- 

Based on the above results, the proposed solution to the reconstruction problem is as follows (Section [3?2]): 

(1) If /C > 0, replace D with 5 =(nf=i(a; - I) • 

(2) Build the matrix H and solve the system _ffx = where x is the vector of unknown coefficients of 2) 
or S according to the previous step. Obtain a solution a and build the differential operator S)* = S)a 
which annihilates g in its entirety (as a distribution in case > or in the usual sense otherwise). 

(3) If AC > 0, recover {^i} and the operator (which annihilates every piece of g) from D* . 

(4) Compute the basis for TV-gt and solve the system Ca = m. 

Conditions for the solvability of the above systems are discussed and few initial results in this direction are 
obtained in Section 13.31 below. Results of numerical simulations, presented in Section 21 suggest that the 
method is practically applicable to several different signal models, including piecewise constant functions 
(and piecewise polynomials), rational functions and piecewise sinusoids. 

1.2. Acknov^rledgements. The author wishes to thank Yosef Yomdin for useful discussions and remarks. 
Also, the criticism of the reviewers has been very helpful. 

2. Recurrence relation for moments of piecewise D-finite functions 

2.1. Single continuity interval. We start with the case of a Z?-finite function g ow a. single continuity 
interval [a, b]. The main tools used in the subsequent derivation are the discrete difference calculus and the 
Lagrange identity for a differential operator and its adjoint. Let us briefly introduce these tools. 

Let s : N ^ M be a discrete sequence. The discrete shift operator E is defined by Es(n) = s{n + 1), and the 

dcf 

forward difference operator by A = E — I. We shall express recurrence relations in terms of polynomials in E 

or A. For example, the Fibonacci sequence Fk satisfies E^ i^^ = E i^^ + Fk, so the operator P(E) E^ — E — I 
is an annihilating difference operator for F^. Likewise, the operator A annihilates every constant sequence. 

Lemma 2.1 ([13j). Let p(Ei) be a polynomial in the shift operator E and let gin) be any discrete function. 
Then 

p(E)(6"g(n)) = 6>(feE).9(n) 

3 



Lemma 2.2 ('13'). For any polynomial p{n) of degree k and for i > 1: 

A'^+Xn) = 

Definition 2.3. Let a; £ R and fc G N. The fcth falling factorial of x is 

{x)k '= x{x - I) . . . [x - k + 1) 

The following well-known properties of the falling factorial are immediately derived from Definition 12.3 
Proposition 2.4. Let x G M, fc G N. Then 



(1) {x)k is a polynomial in x of degree k (thus it is also called the factorial polynomial^. 

(2) If X eN and n > k, then {n)k = j^^,- 

(3) Ifx^neEU {0} and n < k, then {n)k = 0. 

The formal adjoint of a differential operator D = '^f^QPj{x)d^ is given by 

N 

(2.1) ^^w^'Y.i-^y^'iM^y} 

The operator and its adjoint are connected by the Lagrange identity f [19j): for every u,v G 

(2.2) uS)(u) -mD*(u) = -^P2,(u,w) 

dx 

where Px) (""j v) is the bilinear concomitant - a homogeneous bilinear form which may be written explicitly 
as ([B p.211]): 

Ps(w, v)^u {p^v - d{p2v) + ■■■ + (-i)^-ia^-i(pjv«)} 
+ u' {p^v - d{p,v) + ■■■ + i-if-^d^'-^PNy)} 

If (12. 2p is integrated between a and b, Green's formula is obtained: 
(2.4) {I)u,v)-{u,D*v)^[P^{u,v)t 
where the inner product is defined by(it, v) =^ J^^ uvdx. 



Let D and g be arbitrary. Consider the "differential moments" associated with D: 
By ((2^ . we have 



m?(5) =™fc(S)ff) 



mfe(S)g) = - (.9, S) *(.:'=)) + [P^{g,x'')]^ 

Let us define the following two sequences, indexed by k: 

Mfc = Mfc(D,g) -'(5,S*(a:")) 

The sequence of the differential moments is therefore the sum of the two sequences: 
(2.5) mfe(Dg) =Aifc(D,5)+£fc(D,g) 



^ The Pochhammer symbol is also used in the theory of special functions. There it usually represents the rising factorial 
{x)„ = X ■ (x + 1) ■ . . . ■ (x + n - 1). 
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Using and we have: 



N 



(2.6) 



where 
(2.7) 



j=o ^ 

kj N 

^^^aij{-iy {i + k)jmi+k-]{9) = 6s(fc,E)mfc 
1=0 j=o 



e = e£,(fc,E)'ii'^a,.,n(*'^")(fc,E) 



U^''^\k, E) = (-l)J■(^ + k)j E'-^ 



Now let us return to our main problem. Recall that g is Z?-finite, so let T) g — 0. (12. 5p combined with 
gives 6mfc +6^ = 0. As we demonstrate below, there exists a discrete difference operator £ = 8(E) such 
that f Efe = 0. Multiplying the last equation by this £ from the left (multiplication being composition of 
difference operators) gives us the desired recurrence relation: £ ■ Qnik — 0. 

The sequence is related to the behavior of g at the endpoints of the interval [a,b]. The following lemma 
unravels its structure. 

Lemma 2.5. Let D be of degree N as in (|1.2p . Then there exist polynomials qa{k) and qb{k) of degree at 
most N — 1 such that 



(2.8) 



efc(S),5) = b'^qbik) ~ a\a{k), fc = 0, 1, ■ 



Proof. Write D = X^jLo ®i where Dj = pj{x)dK Denote Skj =^ efe(Dj, g). By (|2.3p we have Sk = SjLo 
where 

sk,, = PvAd,^") = [g^'-'^^'p, ~ g'^'-^'^dix^p,) + ■■■ + {-ly-^gd^-^x'^p,) 

Use Leibniz rule: 



(=0 



axx'=p,(x)) (; (x^-)«p(^-')(x) (;)pf')(x)(fc),x'=^ 



i=0 



where rij^i{x) — x ''{fjPj {x) is a rational function. Now 



Sk.j = <x 



'^Y.{~irg^^-'~^\x)Y^{k)ir,.,,{x) 



4=0 

uk „_ n„\ „fc. 



1=0 



1=0 



i=0 1=0 



b^qb^jik) - a'^Qajik) 



where Sijj{x) ^ {-lYg'-^ ^ *)(a;)rij,;(x) and ^^(fc) = X^Lo I]i=o('=)'*ij,'(") Q;e{a,6}. These are 
polynomials in k of degree at most j — 1 (see Definition 12. 3p . Now 

N N N 

j=0 ]=0 j=o 

Take ^^(fc) =^ J2 jLo la.. and (jf,(fc) =^ Ejlo ^i'^*^^ deg(?aj,9bj < j, then we have deg qa,qb < N 

and this completes the proof. ■ 



As a side remark, we have the following simple condition for the sequence {ek} to be a nonzero sequence. 
Theorem 2.6. Assume g = and pn{x) ^ on [a,b\. Then ek{'^,g) = if and only if g = Q. 
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Proof. From the proof of Lcmma l2.5l we have for a G {a, b}: 

(2.9) g„(fc) = ^5:(-l)V^-i-)(a)5:(fc),a-'r 

j=o i=a 1=0 ^ ^ 

• In one direction, we have g{a) — ■ ■ ■ — ^'^^^■'(a) = for a = a,b. By direct substitution we obtain 
qa{k)=0. 

• To prove the other direction, assume b^qb{k) = a^qa{k) — C . Consider the foUowing cases: 

(1) C i.e. a, 6 ^ and qa{k), qb{k) ^ 0. Then (|) ^ = for aU fc e N. But this is impos- 
sible because the left hand side is of exponential growth and the right-hand side is of at most 
polynomial growth. 

(2) C — 0. Then at least one of qa{k),qb{k) must be identically zero. So let qa{k) = and a ^ 0. 
The coefficient of the highest order term of qa corresponds to j — N,i^l = N— 1 and equals 
to 

(-l)^-ig(a)a-(^-iW(a) 

This is zero only if g{a) — 0. So we can lower the limit of the second summation in (|2.9p 
to n — 2. Repeating this argument with j = N,i = I ^ N — 2 shows that g'{a) ~ 0. Finally 
we obtain that g{a) = ■ . . = g^'^~^\a) — and by the uniqueness theorem for linear ODEs we 
conclude g = 0. ■ 

An immediate corollary is that generically we have £ y^l (generically here means g 0). 

With the structure of the sequence {£&} at hand, we can now explicitly construct an annihilating difference 
operator for it. 

Definition 2.7. Given j, a and b, let £ij,{E) (E -alp (E -61)^'. 

Theorem 2.8. Let — 'Yl!j=oPji^)'^^ ■ Then the difference operator £ — £^fj annihilates the sequence 
eki^,g). 



Proof. By Lemma [^31 Ek ~ b^qhifi) — a'^qa{k) where degga, qb < N — 1. The factors (E— al) and (E— 61) 
commute, so it would be sufficient to show that (E — a I)^{a'^gci(fc)} = 0. 

Let P(E) = (E-al)^, then P(aE) = A^. By Lemma [O P{E){a''qa{k)} = a''+^ qa{k). Since 
deg qa < N - I, hy Lemma we have qa{k) =0. ■ 

Remark 2.1. Theorem 12.81 defines a connection between the moments of the two functions g and Tig: 
£ mk{D g) = £Qnik{g). This implies a connection between their moment generating functions as well (see 
further Section [5TT] and also [21]). 

Remark 2.2. The polynomial nature of the coefficients Pj{x) of S has not been used in the proof of Theorem 
12.81 Therefore it is true for every linear operator J) with sufficiently smooth coefficients pj (x) . 

Now we have all the necessary information in order to prove the main result of this section. 

Theorem 2.9. Assume that Tig = Q on [a,b]. Then the sequence {mk{g)} satisfies the following recurrence 
relation: 

(2.10) 6 fUk • e^) mfe = ( (E -aI)^(E -61)^ ^ + fc). E*"' ) = 

\ j=0 i=0 ' 



Proof. We have mk{T) g) = 0. The proof is completed using (|2.5p . (|2.6p and Theorem 12.81 
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Remark 2.3. With respect to {nik}, the length of the recurrence relation (|2.10p is at most 3A^ + ma.xkj + 1. 
Its coefficients are hnear in j and polynomial in k. 

Remark 2.4. The recurrence relation (|2.10p is not trivial, i.e. 6/0. To see this, let J) g = 0. Now if g ^ 
then by Theorem l2.6l £t. ^ and thus /Xfc = 9 ^ 0. It follows that Q cannot be the identical zero operator, 
and therefore & = £ Q ^ 0. 

We conclude this section with some examples which demonstrate the usefulness of Theorem 12.91 

Example 2.1. In [B], the authors provide explicit recurrence relations satisfied by the moments of the 
powers of the modified Bessel function f{x) — Kq(x). The method used to obtain these recurrences is 
integration by parts. However, an additional condition is imposed - namely, it is required that the integrals 
/p x'^(/(a;)")^^'da; converge and the limits of the integrands coincide at the endpoints of T. In our setting, 
this is equivalent to putting £ = 1. 

The function g — Kq{x) is annihilated by the operator 'D — x^d^ + 3xd^ + (1 ^ 4a;^)9 — 4x1 ([Bi Example 
2]). The only nonzero coefficients are therefore 





i = 


i = 1 


i = 2 


j = 




-4 




j = l 


1 




-4 


J = 2 




3 




J = 3 






1 



By (|2T0ll we have 

—Amk+i — knik-i + 4:(k + 2)mk+i + 3fc(fc + l)mk-i — (fc + 2)(fc + l)kmk-i — 

which is just 4(fc + l)mk+i — k'^mk-i- This result agrees with 6, Example 3]. H 

Example 2.2 (Example ll . ll continued) . / = ae^^ is annihilated by the operator D = d — (31. Thus ao,o — —P 
and aifi — 1. By (|2.10p the sequence {rrifc} on [0, 1] satisfies 

6 = E(E-I)(-/3I-fcE-i)mfc = -/3mfc+2 + {/3 - {k + 2))mk+i + {k + l)mk = 

This recurrence relation with polynomial coefficients may be solved explicitly using computer algebra tools 
(see [34] for an overview of the existing algorithms). Using the Maxima computer algebra system (p8j), the 
following explicit formula was obtained: 

Cl = mo 

C2 = (3mi + niQ 

This example is further continued in Example 13.11 I 

2.2. Piecewise case. Until now we have been considering a function g which satisfies S 5 = on a single 
interval [a, b]. In particular, we have seen that the sequence of the moments of g satisfies a linear recurrence 
relation whose coefficients linearly depend on the coefficients of 3D. Now we are going to consider the piecewise 
case (depicted in Figure [1]) where g is assumed to consist of several "pieces" go, gi, . . . , gjc- On each continuity 
interval A„ = [CruCn+i] the n-th piece of g satisfies Dngn — where J)„ — J2j'=o^i=oi'^iJ,'n^^)^'^ ■ Denote 
mk,n J^"^^ x'' gn{x)dx. Our goal is to find a recurrence relation satisfied by the sequence nik = X]n=o ^k,n- 

As we shall see below, this recurrence relation is particularly easy to find explicitly in the case g G AJ,, i.e. 
Dn = S for n — 0, . . . , K.. We shall also discuss the general case where 'Dm 2D„ briefly. 
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Figure 1. Piecewise Z?-fimte function 



2.2.1. T/ie same operator on each interval. We present two methods for computing the desired recurrence 
relation. Then we show that these methods indeed produce the same result. 

Theorem 2.10 (Method 1). Let g e AJ, with Dm = D for < m < K.. Then {mk{g)} satisfies 



(2.11) 



K.+1 



[](E-C„I)^es(fc,E) TOfe = 



Proof. By Theorem 12.91 for every n = 0, . . . , /C the sequence {mfc.n} satisfies 

(2.12) (^6.,«„+i(E) • es,(fc,E))mfe,„ = 

p. lip follows by using the fact that the linear factors (E —^i I) and (E — I) commute. 



It turns out that the formula (|2.1ip may be obtained by considering the piecewise function g G AJ, being an- 
nihilated as a distribution on [a, b] by some operator S) = X)(S, {Cn})- We now derive the explicit expression 
for S. 

The general theory of distributions (generalized functions) may be found in [16] . By a test function we shall 
mean any / 6 C^~^{[a, b]) (in fact, for our purposes it is sufficient to consider just the moments / — x^). 

We shall identify the discontinuous function g with the following distribution (here g_i =0 by definition): 

K 

(2.13a) g{x) ^ go + ^'g^{x)n{x - £,n) 



(2.13b) 
(2.13c) 



~ dcf 

Qn — gn ~ gn-1 



nix] 



dcf 



X < 

1 X > 



The functions gn{x) belong to the solution space of D on A„, and thus gn{x) G ^{An). Since Pn{x) 7^ 
on [a, b], we even have g^{x) G C^~^{[a, b]). 

The derivative of H is the Dirac 6. The distribution S{x) has the following properties: 

{S{x-t),f{x))=f{t) 

{s<^^\x-t)j{x)) = {-iyf^^\t) 



The second equality is valid provided /(x) is j-times differentiable at x — t. 



Now let g{x) — c ■ S{x — ^) for c, ^ e M. This distribution is "annihilated" by the operator D = (a; — ^) I 
in the sense that (2) g, /) = (c • (x — ^)5{x — £,),/) =0 for every test function /. By linearity, the operator 
S) — Y[f=ii^ ~ ^i) I annihilates every g of the form g — Ylf=i Ci5{x — ^i). 

Lemma 2.11. Let g — X)J=o^ Cij{x)S^^^ {x — ^i) such that Cijix) is j -times differentiable at ^i. Then g 

is is annihilated by J) = p{x) I where p{x) = Y\f^i{x — ■ 

Proof. For every test function f{x) we have 

K S-l 

(/, Dg) = {p{x)f{x), 9)^Y.T. -^^'^ - 

i=l j=0 

The function r{x) — Cij{x)p{x)f{x) has a zero of order 5 at a; = Therefore, all its derivatives up to j 
vanish at £,i and so (/, D g) — Q. ■ 

Theorem 2.12. Assume g S AJ, as in (|2.13p . with X) of degree N annihilating every piece of g. Then the 
entire g is annihilated as a distribution by the operator D = JliLi(^ ~ ^i)^ ®- 

Proof. By (|2.13p and the fact that D = we obtain 

N r K. 



IC N j 



1=1 j=0 k=0 ^ ^ 

1=1 j=0 fc=l ^ ^ i=l j"=0 

K Af-1 K 



i=l k=Q 1=1 

= EE'^^^(^)'^*'^(^-co 

i=l fc=0 

where hik{x) is fc-times differentiable at Then apply Lemma 12.111 ■ 

Theorem 12.121 will later serve as the basis for recovering the locations of the discontinuities of G AJ, . The 
factor Y\{x — ^i)^ effectively "encodes" the positions of the jumps into the operator itself. The "decoding" 
will then simply be to find these "extra" roots, once the "enlarged" operator is reconstructed. 

The second form of the recurrence relation now immediately follows from Theorem 12.121 

Theorem 2.13 (Method II). The sequence of the moments of g £ AJ, satisfies the recurrence relation 

(2.14) (f^,(E).eg(A:,E))mfe = 

where!) ^l\n=ii^-^n)^ ^■ 

Proof. Theorem 12.91 combined with Theorem 12 .121 ■ 
Lemma 2.14 (Equivalence of Methods I and II). Let q{x) be an arbitrary polynomial. Then 

0,(.)s(A;,E) = 9(E)es(A;,E) 



Equivalence of (|2.1ip and (|2.14p then follows from Lemma [2. 141 by putting q{x) = 11^=1 ^ Ci)^- 
To prove Lemma 12.141 we need the following result. 

Proposition 2.15. Let p{x) = ^f^oPiX^, q{x) = X]f=o ^i-^' andr{x) = p{x)q{x) — ^^"=0 ^i^^ ■ LetE be the 
shift operator in k and let j be fixed. Then 

Q+/3 ^13 _ 

^ r,(* + fc)^E''-^ = p{E) ^ + k)jE^-^ 

i=0 i=0 

Proof. We extend the sequences of the coefficients {pi} and {qi} by zeros as necessary. By the rule of 
polynomial multiplication, = J2'j=oPj9i-r 

Q+/3 Q+/3 a a Q+/3 

(■.-.-0 = T.P^ E^" E + fc)7E'-^ = P(E) E ^^(^ + fc)7E'"^ ■ 

Proof of Lemma l2l4\ Recall that 6© = X^ij aijn^^^^) (fc, E) where H^*'^) (A:, E) = (-l)J(i + fc)j E'"^ 2) = pj(a 
and Pj{x) = 'J2iLo Now let Pj{x) = q{x)pj{x) — Y^ai^jX^ , then 

AT fej+dcgij 

e,(.)s(fc,E) = E(-i)' E «m(* + %e^"'- 

j=0 1=0 
AT kj 

(ProvosiUon rSTSi = E (~ ^'^ ^^^^ E "^'-J'^* ^ ^'"'^ 
i=0 i=0 

= g(E)e»(fc,E) ■ 

2.2.2. Different operators on each interval. Recall that we want to find a recurrence relation for the sequence 
m/c = ^ rrik^n where each subsequence is annihilated by the difference operator 

®« = ^l!:«„+rQs„(fc,E) 

There exist at least two approaches, both of which involve techniques from the theory of non-commutative 
polynomials - the so called Ore polynomial rings (see [Hj). Both differential and difference operators with 
polynomial coefficients are members of the appropriate Ore algebra. The least common left multiple (LCLM) 
of two polynomials p, q is the unique polynomial r of minimal degree such that both p and q are right-hand 
factors of r, i.e. r = p'p = q'q for some polynomials p', q' . The LCLM may be explicitly found by the 
non-commutative version of the polynomial division algorithm. The complete theory may be found in [26j . 

Approach I Given the operators S„ such that S„ TOfc,n = 0, the operator &^ which annihilates the sum 
™fc,ri is given by the least common left multiple of {S„}. 
Approach II Given the operators £)„ which annihilate the pieces g„ separately, the operator which 
annihilates every piece simultaneously is the least common left multiple of {£>„}. Then the 
annihilating operator for {nik} is = £^f, Oj,t(A:,E). 

We are not aware of any general procedure by which the coefhcients of 6^ or &^ may be related to the 
coefRcients of each S)„ in some tractable manner, unless the operators S)„ commute. 

Example 2.3 (Piecewise sinusoids). Let g consist of two sinusoid pieces: gi{x) = ci sin((jjia; + ^i) and 
92{x) = C2 sin (0^22; 4- 02) with break point ^. The annihilating operators for gi and 172 are, respectively, 
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Di = d'^ + ujfl and D2 = d'^ + u)\ I. The operator D ' = (5^ + ^^i I) • (5^ + lJ\ I) annihilates both pieces si- 
muhaneously, therefore [x — annihilates the entire g. H 



3. Moment inversion 

We now present our method for moment inversion for piecewise D-finite functions. Recall that our purpose 
is to reconstruct the parameters 



(3.1) 

from the input 
(3.2) 



^ - {K,},{e.},K„}} 



X =^ {{mfc},7V, {fcj},/C,a, 6} 



First we establish explicit connections between V and X - the "forward mapping" : P — > I. Then we 
derive the inverse mapping M ~ M.^^ and provide simple conditions for the solvability of the resulting 
inverse systems. 

3.1. Forward equations. Recall the polynomials n(''^^(A;, E) which were defined in (|2.7p during the deriva- 
tion of the recurrence relation (|2.10p . Given a multiindex (i, j) we define the "shifted" moment sequence 

(3.3) .f'^=^^(<,(E).n('.^)(fc,E))m, 
For each j = 0, . . . , let hj(z) be the formal power series 



fc=0 



Finally, for any g(x) let 



00 



k=0 



Ig{z) is called the moment- generating function of g. 

Theorem 3.1. Let g G A'^ he annihilated by S (either in the usual sense if IC = or as a distribution if 
IC > 0, in which case D = Y\f^i{x — ^i)^J)^ with J)^ annihilating every piece) , where D = X^^LoPjl^)^"' 



The 



(A) The vector a — (aij) satisfies a linear homogeneous system 



(3.4) 



Ha 



,(0,0) „,(1,0) 



(0,0 (1,0 



(kN,N) 


1 



M 



M 



^(fe„,7V) 
M 



( "0,0 \ 



oi 



\afejv,JV/ 



for all M e N. 

(B) v^*'"''' ~ rni^k {Sjix)) for all < j < N , < i < kj and fc G N (the moments are taken in [a, b]). Con- 
sequently, hj{z) is the moment generating function of q^-{x). 
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( C) The functions {1, ho{z), . . . /lAr(z)} are polynomially dependent: 



N 



where Pj{z) — z™^^^^ pj{z ^) and q{z) is a polynomial with degq < maxfcj. 
Proof o/Hl By Theorem [^J^ and f^ we have 

AT fej 

EE«^^^P = 0' ^ = 0,1,... 

This is exactly 

Proposition 3.2. Let p{x) be a polynomial in x. Then for every f{x) 

ink{p{x)f{x)) ^ p(E)mk{f{x)) 

Proof. Let p{x) = x^ . Then 

mk{x'' f{x)) ^ I x^x''f{x)dx^ I x^+''f{x)dx^mk+r{f{x))^Wmk{f{x)) 

J a J a 

The proof for an arbitrary polynomial p{x) follows by linearity. 

Proof o/O Fix j < N, i < maxfcj and define — x^d^ . By p.5p and (|2.6p we have 
(3.5) mfc(£>y g) = 63,^. (fc, E)mfc(5) + £^(2),^, 5) 

By dlJl) we have: es^^(fc,E) = H^^'^) (fc, E). By Theorem[2l]we have f ^b(E)efe(Dy , .g) = 0. Finally, 



<,(E) • n(^'^") (fc, E))mfc(g) = f ^,(E) 6^,, (fc, E)™,. 



|33J = fl(E)mfc(D,,- g) - £:^,(E)£fc(S).„ 5) - <b(E)™fc(S)., g) 
(Pr-opo5»ton[UJ = mk{£{x)x'd^g) = TO,+fc(g^(a;)) 

Proof 0/0 Let fc* = maxfcj. We have Pj{z) = YliLo o-i^jz'^ ~*. Denote the power series 

N 00 
j=0 fc=0 

An immediate consequence of[B]is that u^! = w^* "''' for i' + fc' = i + fc and all j. Then for all m > k* 







j=0 1=0 j=0 i=0 

So g(z) is a polynomial of degree at most k* — 1. This completes the proof of Theorem 13. II ■ 

Remark 3.1. H has the structure of Hankel- striped matrix H = [Vq . . . V^] where each "stripe" is a Hankel 
matrix 



(3.6) 



■M (ij) 
(O-j) (i.j) 



■ M Af 



-'o 



M 



Hankel-striped matrices appear as central objects in contexts such as Hermite-Pade approximation (the 
standard Pade approximation being its special case), minimal realization problem in control theory and 
Reed-Solomon codes ( 23, 20, 7 ). In fact, the system of polynomials {Pj{z)} is called the Pade-Hermite form 
foT<P^{l,ha{z),...,hNiz)}. 

Remark 3.2. Moment-generating functions are a powerful tool for the investigation of the properties of the 
sequence nik- For instance, the asymptotic behavior of the general term may be derived from the analytic 
properties oi Ig{z) ([T5]). 



Now suppose the operator S) annihilating every piece gn of g £ AJ,, as well as the jump points are 
known. Let {wijf^^ be a basis for the space Afs- Then gn{x) = '^^i (^^^{{x) . Applying the moment 
transform to both sides of the last equation and summing over n = 0, . . . ,IC gives 

Proposition 3.3. Denote c^^. = /|^"+' x''u^{x) /or n = 0, . . . , IC. Then VM e N; 



(3.7) 



V' 



l.M 



N,M 



-N 



A 



N,M/ 



( "1,0 \ 



\aN,K/ 



nil 



L 



Input: {mfcj^l^o' ^>{^j}> ^ 




Yes 



Reconstruct £> 
Solve Ha = (lOl) 




Yes 



No 



7 



1 



,s)-nr=i(^-6)^s 



J 



Recover jump points 



Recover particular solution(s) 
Find basis for TVVit . Solve Ca — m (13.71 



Figure 2. The reconstruction algorithm 
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3.2. The inversion algorithm. The reconstruction algorithm which is based on the results of the previous 
section is depicted schematically in Figure [2l The solvability of the corresponding systems is discussed in 
the next section. Note the following: 

(a) At the initial stage, the "encoding" of the (yet unknown) jump points takes place. In practice, this means 

the "enlargement" of J) to J) =^ niLi(^ ~ f*)^ ®- The parameters of the problem therefore change as 
follows: N remains the same while kj 4— kj + NJC - see (|2.14p . 

(b) By Theorem 12. 121 {^n} are the K. distinct common roots of the polynomials which are the coefficients 
of D of multiplicity N. The remaining part of the coefficients define the operator J)^ which annihilates 
every piece of g. 

Example 3.1 (Examples 11.11 and 12.21 continued) . g{x) = ae^^ on [0, 1] is annihilated by £) = d — f3l. 
Writing down ([23) with M = yields 

[m2 - mi -{2mi - mo)] 

which has the solution 

^ _ 2mi - uiQ 
mi — m2 

The constant a is then recovered by 

mo{g) f3mo 
d — = 

mo(e'^^) e'^ — 1 

Note that this solution requires the first 3 moments instead of 2 as in (|1.4p . ■ 
Additional examples of complete inversion procedures are elaborated in Appendix [Xj 



= 



3.3. Solvability of inverse equations. The constants M and M determine the minimal size of the corre- 
sponding linear systems (|3.4p and (|3.7p in order for all the solutions of these systems to be also solutions of 
the original problem. 

Theorem 3.4. //b <^Mh, then 

mk{E{x){^^,g){x))^^ k = 0,l,...,M 



Proof. Denote b ~ {bij) and let A; > 0. By Theorem 13.11 Part IB] and Proposition 13.21 the product of the 
k + 1-st row of H with b is 

= EE^^^^i^''^ = ^X^6„m,+,(flf(a;)) ^ ^ m,(6,,:rX^(x)) = m, E(E^»^^>'(^) 

j=0 i=0 j=0 i=0 j=Q i=0 \j=0 ^i=0 ' ) 

= mfe E{x){T)t,g^{x) 



If it is possible to estimate how many moments of _F = E{x){^^a g)(x) should vanish in order to guarantee 
the identical vanishing of F and therefore also of IDb g for all possible differential operators of the prescribed 
complexity (i.e. of given order and given degrees of its coefficients), then M may be taken to be this number. 
Then, every solution of (|3.4p will correspond to some annihilating operator. For any specific g and such a 
finite number exists, because any nonzero piecewise-continuous integrable function has at least some nonzero 
moments. However, this number may be arbitrarily large as shown by the next example. 
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Example 3.2 (Legendre orthogonal polynomials). It is known that every family of orthogonal polynomials 
satisfies a differential equation of order = 2 of the following type: 

(3.8) ^ ^ q{x)d'^ + p{x)d + Xnl 

where q^p are fixed polynomials with degq < 2, degp < I and A„ is a scalar which is different for each 
member of the family. 

Consider {Ln{x)} - the family of Legendre orthogonal polynomials. The interval of orthogonality is [a, b] — [—1, 1] 
and £ = (E^ — I)^. i„ is annihilated by £>„ = (1 — x^)d'^ — 2xd + A„ I. Furthermore, (L„(x), x'^) = for 
k < n — 1. The "reconstruction problem" for Ln{x) (assume it is normalized) is to find the constant A„ from 
the moments (it is well-known that A„ = n(n + 1)). 

Take an arbitrary vector b = [6oo bu 602 ^22]^ so S)b = bool+buxd + (&02 + b22X^)d^. The function 
£{x) Db Ln is a polynomial of degree n + 4, so it is uniquely determined by its n + 5 moments (t27j). By 
Theorem 13.41 this would be the minimal size of the system p.4p . The entries of the fc-th row of H are 

"fc°'°^ = riik+i - 2mfe+2 + mk i = 0, j = 

v'^^''^^ = -((/;; + 5)mfe+4 - 2{k + S)mk+2 + (k + l)mk) i = 1, j = 1 

4°'^^ = {k + ^){k + 3)mk+2-2{k + 2){k + l)mk + k{k-l)mk-2 i = 0,j = 2 

^ (fc + 6)(fc + 5)mfe+4-2(fc + 4)(/s + 3)mfc+2 + (fc + 2)(/c+I)mfe i^2,j = 2 

The first n — 5 rows are identically zero because the maximum moment of L„ involved is m„_i. We are 
looking for a solution of the form b = [A„ —2 1 — l] • The {n — 4)-th row is 

[m„ -(n+l)m„ (n + 2)(n + l)m„] 

and that is just enough to reconstruct the operator: m„(A„ + 2(n + I) — (n + 2)(n + I)) = and thus 
A,j = n{n + 1), as expected. The reason for the original estimate n + 5 not being sharp is that some of 
the coefficients of J) (in fact, all but one) were known a priori. ■ 

Similar argument may be used in order to estimate M . Straightforward computation leads to 
Proposition 3.5. If a = {ain) is a solution of p.7p . then 

(3.9) mfc(©„(.x)) -0, A: = 0,I,...,M 

where 

K N 

(3.10) &ol[x)'^= g{x) -'^'^ainUi[x) ■ 

Every function ©a(a;) is a (piecewise) solution of 33 / = 0. By Theorem 12 . 131 the moments of <&a{x) satisfy a 
linear recurrence relation (|2.14p . Therefore, the maximal number of moments of ^oi{x) ^ which are allowed 
to vanish is explicitly bounded by (|2.14p . For instance, M may be taken as the length of the recurrence plus 
the value of the largest positive integer zero of its leading coefficient. 

Example 3.3 (Piecewise-constant functions on [0, 1]). Each piece of 5 is a constant gi{x) = Ci. The operator 
J) — ni=i(^ ~ Ci)^ annihilates entire g and £(E) — E(E — I). 

• First we determine the minimal size of the system p.4p . For every polynomial q{x) of degree /C, the 
moments of the function fq = {£{x)q{x)d) g are 

K 

mkifq) = J2 * = ~ ^)«(?*) 

1=1 

15 



Assume that the ^,;'s are pairwise distinct and do not coincide with the endpoints. We claim that 
M = fC — 1 is sufficient. Indeed, assume that mk{fq) = for fc = 0,l,...,/C — 1. Then 







1 


1 


1 






a 


6 • 






X ^ 




CI • 


■ 






>1 


'52 





92 



The matrix X is a Vandermonde matrix and it is nonsingular because ^ ^ for i ^ j. It foUows that 
q = and therefore q{S,i) = for all i. Therefore every solution of (|3.4p is a multiple of J|,^-^(a; — a). 
Now we determine the minimal size of p.7|) . The space A/© is spanned by piecewise-constant 
functions with the jump points ^i. Let a be a solution to p.7p and let the function <3cy.{x) be as 
in p.lOp . This 0Q.(a^) is again a piecewise-constant function with /C jump points ^i. Therefore, 
the moments of i&a{x) satisfy the recurrence relation (|2.14p . This recurrence has nonzero leading 
coefficient and its length is ^ + 2. Therefore, vanishing of the first IC + 1 moments of &a{x) implies 
= and consequently M = IC will be sufficient. The constants Ci are precisely the solution of 
(I3J1). ■ 



4. Stability of inversion in a noisy setting 



The presented inversion scheme assumes that the unknown signal g G AJ, is clean from noise and that 
the moment sequence {nik} is computed with infinite precision. In many applications this assumption is 
unrealistic. Therefore, it is practically important to analyze the sensitivity of the inversion to the noise in 
the data, both theoretically and numerically. In Sect ion |4T] below, we give a theoretical stability estimate for 
one special case of linear combinations of Dirac (5- functions. Then we present results of numerical simulations 
for several test cases (Section 14. 2p . 



4.1. Theoretical stability analysis. The general case of an arbitrary piecewise _D-finite function appears 
to be difficult to analyze directly. Here we provide stability estimates for the reconstruction of the model 

(4.1) .g(a;)==^a,(5(a;-C.) 



We argue that it is crucial to understand the behaviour of the reconstruction in this special case. Consider 
a generic g £ . Then D g is a linear combination of (5- functions and their derivatives (see the proof of 
Theorem I2.12p . Thus the model (|4.ip may be considered a "prototype" which captures the discontinuous 
nature of g. 

One can prove the following 

Theorem 4.1. Let g be given by (j4.ip . Assume that the moments mk{g) are known with accuracy e. There 
exists a constant Ci such that the parameters may in principle be recovered with the following accuracy: 

|Aa,| < Cie 

where Ci depends only on the geometry of ^i, . . . , ^k:- 
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Outline of the proof: Write the Jacobian matrix of A4 explicitly and use the inverse function theorem to get 
the Jacobian of J\f. This matrix may be factorized as a product of a diagonal matrix diagjai, . . . , a/c} and 
a Vandermonde matrix V on the grid ^i, . . . ,£,ic- Then Ci may be taken to be the norm of V~^. M 



The above result can be generalized to the case of the model 

i=l j=0 

We plan to present these results in detail separately. 

The estimate of Theorem 14.11 reflects only the problem conditioning (i.e. the sensitivity of the map Af defined 
in Section [3] with respect to perturbations). Moreover, this estimate is sharp in the worst-case scenario, since 
it is based on directly evaluating the norm of the Jacobian. 

On the other hand, an accurate analysis of the algorithm itself should involve estimates for the condition 
numbers of the matrices H and C (see ()3.4p and (|3.7p . respectively) as well as the sensitivity of the root 
finding step (Figure [2). Let us briefly discuss these. 

(1) By Theorem 13. II and Remark |3. 11 solving p.4p is equivalent to calculating a Hermite-Pade approx- 
imant to $ = {l,ho{z), . . . jh^^z)}. Estimates for the condition number of H in terms of <i> are 
available ([8]). These estimates may hopefully be understood in terms of the differential operator D 
and the function g itself, using the connection provided by Theorem 13. II 

(2) It is known (e.g. [31]) that if the coefficients of a (monic) polynomial are known up to precision e, 
then the error in recovering a root of multiplicity m may be as large as e " . 

(3) There is some freedom in the choice of the basis for the null space of D (Proposition 13. 3p . 
It should be investigated how this choice affects the condition number of C. 

Considering the above, we hope that a stable reconstruction is possible at least for signals with relatively 
few and sufficiently separated points of discontinuity, as well as some mild conditions on the other model 
parameters. 



4.2. Numerical results. In order to further justify the hope expressed in the previous section, we have per- 
formed several numerical simulations using a straightforward implementation of the reconstruction algorithm 
from Section [3] (the implementation details are provided in Appendix . 

We have chosen 3 different models for the unknown signal: a rational function (no discontinuities), a piece- 
wise sinusoid, and a piecewise-polynomial function. In every simulation, the signal was first sampled on a 
sufficiently dense grid, then a white Gaussian noise of specified Signal-to-Noise Ratio (SNR) was added to 
the samples, and then the moments of the noised signal were calculated by trapezoidal numerical integration. 
To measure the success of the reconstruction, we calculated both the Mean Squared Error (MSE) as well as 
the error in the relevant model parameters. The results are presented in Figures [3J U] and [S] Note that the 
SNR is measured in decibels (dB). 

Several test signals were successfully recovered after a moderate amount of noise was added. Furthermore, 
increase in SNR generally leads to improvement in accuracy. Nevertheless, some degree of instability is 
present, which is evident from the peaks in Figures [3(d)[ |4(by[ |5(d) 



While for some models, the increase in complexity leads to severe performance degradation, for other models 



it is not the case. Compare the growth in the degree of the rational function on one hand (Figure 3(c) I, 
and the increase in the number of jump points for the piecewise-constant reconstruction on the other hand 
(Figure 5(c) I . 
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Reconstruction of rational signal f=p/q 
degree=[2,4],SNR=25.0251, MSE=0.0710057 



Reconstruction of rational signal f=p/q 
[3,1],SNR=35.1012, I^SE-0.000602905 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



(a) 



(b) 



Reconstruction of rational signal f=p/q 
SNR=40 



Reconstruction of rational signal f=p/q 
degree=[2,3] 





(c) 



(d) 



Figure 3. Reconstruction of rational signals f—p/q- (a) The signal with 
degp = 2, deg q = A corrupted by noise with SNR=25 dB. Reconstruction MSE is 0.07. (b) The sig- 
nal with degp = 3, degg = 1 corrupted by noise with SNR=35 dB. Reconstruction MSE is 0.0006. 
The reconstructed signal is visually indistinguishable from the original, (c) Dependence of the MSE 



on the degree of the denominator with SNR=40 dB. (d) Dependence of the MSE on the SNR with 
degp = 2,degg = 3. 



5. Discussion 



The piecewise D-finite moment inversion problem appears to be far from completely solved. In the theoretical 
direction, the solvability conditions need to be further refined. In particular, we hope that minimality results 
("what is the minimal number of measurements which is sufficient to recover the model uniquely?") in the 
spirit of Examples l3.2l and l3.3l concerning the reconstruction of additional classes of functions may be obtained 
using the methods presented in this paper (Section 13. 3p . The importance of this question is discussed e.g. 
in [27| , where estimates on the finite moment determinacy of piecewise-algebraic functions are given. In this 
context, the role of the moment generating function is not yet fully understood. Another open question is 
the analysis of the case in which the function satisfies different operators on every continuity interval. 
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Reconstruction of piecewise sinusoidal signal, pieces=4 
SNR=24.986, MSE=Q.0675575 



Piecewise sinusoid reconstruction, pieces=2 
[a,b]-[-1,11, Jump-0.45, Freq.=3 





(b) 



Figure 4. Reconstruction of piecewise sinusoids. |(a)| A sinusoid consisting of 4 pieces, 
corrupted by noise with SNR=25 dB. Reconstruction MSE=0.068. [(b)] Dependence of the MSE 



on the SNR in the recontruction of a 2-piece sinusoid. Also plotted are the relative errors in the 
estimated location of the jump point and the estimated frequency. 



In the numerical direction, the results of Section [4] suggest that a "naive" implementation of the presented 
algorithm is relatively accurate for simple enough signals corrupted by low noise levels. We believe that 
attempts to improve the robustness of the algorithm should proceed in at least the following directions: 



(1) A similar question regarding stability of reconstruction of signals with finite rate of innovation is 
addressed in [53]. We propose to investigate the applicability of that method to our reconstruction 
problem. 

(2) The connection of the system p. 41) to Hermite-Pade approximation may hopefully be exploited to 
build some kind of sequence of approximants which converge to the true solution, as more elements 
of the moment sequence are known. In fact, a similar approach (with standard Pade approximants) 
has been used in the "noisy trigonometric moment problem" ([3]). 



Furthermore, we believe that it is important to understand how the various model parameters influence the 
stability of the algorithm. The most general answer in our context would be to give stability estimates in 
terms of the differential operator S and the geometry of the jump points. 

On the other hand, the algebraic moments (|l.ip are known to be a non-optimal choice for measurements 
(see |32j ) due to their strong non-orthogonality. As stated in the Introduction, the moment inversion is a 
"prototype" for some real-world problems such as reconstruction from Fourier measurements. In fact, it is 
known that the Fourier coefficients (as well as coefficients with respect to other orthogonal systems) of many 
"finite-parametric" signals satisfy various relations, and this fact has been utilized in numerous schemes for 
nonlinear Fourier inversion ([5] [9] [121 [HI [l] [4]). The point of view presented in this paper, namely, the 
differential operator approach, can hopefully be generalized to include these types of measurements as well. 
In this regard, we expect that methods of holonomic combinatorics (|35j) may provide useful insights. 
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Reconstructing piecewise polynomial signal, degree=0 
SNR=1 4.9907, MSE=0.0432778 



Reconstructing piecewise polynomial signal, degree^l 
SNR=30.0ia, MSE-0.0674048 
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Piecewise polynomial reconstruction 
degree=0,SNR=50 



Piecewise polynomial reconstruction 
degree=0, pieces=3 



-MSE 
Jump Err. 
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Number of jump; 



Figure 5. Reconstruction of piecewise-polynomials. (a) Piecewise-constant signal with 



5 jumps corrupted by noise with SNR=15 dB. Reconstruction MSE is 0.043. The reconstructed 



signal is visually indistinguishable from the original, (b) Piecewise-linear signal with 3 jumps. 
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Appendix A. The inversion algorithm - implementation details 

We have implemented the inversion algorithm of Section [3] in the JVIATLAB environment For each of 

the three types of signals we built a specific inversion routine. All the three routines have a common base 
as follows: 

• The matrix H is chosen to be square. Solution to the system Ha. is obtained by taking the 
highest coefficient of a to be 1 and then performing standard Gaussian elimination with partial 
pivoting on the reduced system. 

• The root finding step is performed with the pejroot routine of the MULTROOT package ([321 137]). 
which takes into account the multiplicity structure of the polynomial. 

• The calculation of the moments in the matrix C is done via numerical quadrature. 
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The difference between the routines lies in the construction of H and the computation of the basis for the 
nuUspace of S). These are determined by the signal type, as described below. 



(1) Rational functions: A rational function /(x) — is annihilated by the first order operator 

D = {-pq)d+{p'q~pq')l 

Therefore the matrix H consists of two "stripes" Vb and Vi (the structure of the matrices Vi is given 
by ()3.6|) ). If Aegp{x) = r and deg q{x) = s, then the degrees of the coefficients of J) are — r ~V s — 1 
and ki — r + s. Subsequently, Vq is (/cq + ^"i + 2) x (fco + 1) and Vi is (fco + ki +2) x (fci + 1). After 
D is reconstructed (without explicitly recovering p and q), a solution m to D u = is obtained by nu- 
merically solving (using the ode45 solver) the initial value problem {J) u{x) = 0, u(a) — l,x — a . . .b}. 
Then the final solution is obtained as 

/» ^ i^.(.) 
mo 

(2) Piecewise-sinusoids: Let g{x) consist of /C + 1 pieces of the form 

gn{x) = An sm{u;x + tpn) 
The annihilating operator for each piece is = + u!^ I, therefore the entire g is annihilated by 
T) = {{x-^,f-....{x-^^f}id^+LuH) 

So H ^ [Vo V2] , where both Vq and V2 are 2(2/C + 1) x (2/C + 1). From the operator reconstruction 
step we obtain D — pi{x)d^ + ^0(2;) I- The jump points are taken to be the arithmetic means of the 
corresponding roots of pq and pi . Because of the normalization performed in the operator recovery 
step, the frequency w' may be taken to be the square root of the highest coefficient of po. The 
reconstruction is considered to be unsuccessful if at least one root fails to lie inside the interval [a,b]. 
The basis for the nuUspace of D is chosen to be ui — sm{uj'x) and U2 = cos{lu'x). 

(3) Piecewise-polynomials: Let g consist of /C + 1 polynomial pieces gn, where degg^ = d. The annihi- 
lating operator for every piece is 5) = d'^^^, and therefore g is annihilated by 

^ ^ {{x ~ ^1)^+' ■ . . . ■ {x - ^^)^+'}d''+' 

Consequently, H ~ Vd+i is a IC{d + 1) x IC{d + 1) Hankel matrix. Operator reconstruction step 
gives D = p{x)d'^^^, and the jump locations are the roots of p{x), each with multiplicity d + 1. The 
reconstruction is considered to be unsuccessful if at least one root fails to lie inside the interval [a, b]. 
The basis for the nuUspace of D is always chosen to be {1, x, . . . , x"^}. 
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